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We introduce a method for treating soft modes within the analytical framework of the quasihar- 
monic equation of state. The corresponding double-well energy-displacement relation is fitted to a 
functional form that is harmonic in both the low- and high-energy limits. Using density-functional 
calculations and statistical physics, we apply the quasiharmonic methodology to solid periclase 
(MgO). We predict the existence of a B1-B2 phase transition at high pressures and temperatures. 
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I. INTRODUCTION 

The quasiharmonic approximation |l[ provides a 
means of extracting finite-temperature properties of ma- 
terials from static calculations. It assumes the vibra- 
tional properties can be understood in terms of ex- 
citations of non- interacting harmonic normal modes: 
phonons. Lattice dynamics g] can be used to calculate 
phonon energies by evaluating the eigenvalues of the dy- 
namical matrix, which involves second derivatives of the 
crystal energy with respect to atomic displacements. 

The frequencies of these modes depend on the crystal's 
density. Hence they have a temperature-dependence that 
arises simply because of thermal expansion in the mate- 
rial. 

Recent developments in ab initio energy calculations 
have enabled full phonon dispersion curves to be ob- 
tained, leading to a resurgence of interest in the quasihar- 
monic approach. Difficulties arise, however, when the dy- 
namical matrix has negative eigenvalues, indicating that 
the crystallographic structure is not a local minimum of 
energy. Such crystals may be dynamically stabilized: be- 
cause of their large entropy, they may represent a mini- 
mum of free energy at high temperature. Within the con- 
ventional assumption of harmonic phonons, the quasihar- 
monic framework leads to divergent free energies in these 
cases. Such systems have been treated numerically from 
first principles using Monte Carlo methods with effective 
Hamiltonians || or molecular dynamics simulations of 
a reduced set of modes Here, we relax the quasi- 
harmonic assumption of harmonic modes while retaining 
the approximation of non-interacting phonons and show 
how the intrinsic anharmonicity of such modes can be 
included in analytic free energy calculations. 

The mineral periclase (MgO) is of some geological 
significance as one of the supposed constituents of the 
Earth's lower mantle. It is generally believed that along 
the geotherm — the conditions of pressure and tempera- 
ture actually occurring in the mantle — periclase remains 
in a single phase. Under other conditions, however, pre- 
vious calculations have suggested that periclase has two 



phases in its solid state: a sodium chloride-like face- 
centered cubic phase (Bl) and a cesium chloride- like sim- 
ple cubic phase (B2) that is favored at extremely high 
pressures ||. 

Imaginary phonon frequencies are found in the B2 
phase of periclase and have attracted enormous atten- 
tion in the other principal constituent of the lower man- 
tle, magnesium silicate perovskite (MgSiOa) [|6| — p~0|] - 

Equilibrium structures, thermodynamic properties and 
compositions depend on the free energy. Here we 
present a calculation of the free energy of periclase as 
a function of density and temperature. We use the 
pseudopotential-plane wave approach to evaluate total 
energies and the method of finite displacements to eval- 
uate pressure-dependent force constants |ll L including 
effective charges and dielectric constants [12] for the 
longitudinal optic modes. Based on calculated phonon 
frequencies and the quasiharmonic approximation, we 
present a first-principles calculation of the phase diagram 
and thermodynamic equation of state of solid periclase: 
the relationship between pressure, density and tempera- 
ture. 



II. THE QUASIHARMONIC METHOD 

A. Ab initio calculation of specific Helmholtz free 
energies 

The first stage of the calculation is to obtain the spe- 
cific (with respect to mass) Helmholtz free energy of each 
phase as a function of density and temperature. We write 
the free energy as the sum of the frozen-ion interaction 
energy and the free energy due to lattice vibrations. 



1. The frozen-ion energy 

The frozen-ion energy — the interaction energy of the 
crystal with the ions fixed in their equilibrium positions — 
is, by definition, temperature-independent. Hence the 
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free energy is simply equal to the internal energy. In or- 
der to determine the dependence on density, total energy 
densiyty functional calculations are carried out for each 
phase at a range of different lattice parameters. 



2. The lattice thermal Helmholtz free energy 

We calculate the lattice thermal Helmholtz free en- 
ergy of each phase as a function of density and temper- 
ature within the framework of the harmonic approxima- 
tion @. 

For a range of lattice parameters, we evaluate the ma- 
trix of force constants for a supercell (several unit cells) 
of the phase under consideration, subject to periodic 
boundary conditions. We evaluate the forces on the ions 
in a crystal when one ion is displaced slightly from its 
equilibrium position: from such calculations the matrix 
of force constants may be constructed |ll[] . 

We denote the matrix of force constants by (j>, where 
4>l,n,a;m,p,p is the component of force in direction a on 
ion n in unit cell I when ion p in unit cell m is displaced 
infinitcsimally in direction /?, divided by the magnitude 
of the displacement. 

Pairs of density-functional calculations are carried out 
with an ion displaced from equilibrium along one of the 
Cartesian axes by a small amount in first a positive 
and then a negative sense. By averaging the resulting 
Hellmann-Feynman forces on the ions from the first sim- 
ulation with the negative of the forces from the second, 
first-order anharmonic contributions to the force con- 
stants are eliminated. 

The set of rotations under which the crystal structure 
is invariant are identified and the rotation matrices, to- 
gether with the mappings between the ions under the 
symmetry operations, are evaluated. For a given pair 
of ions (l,n) and (m,p), the matrix of force constants 
4>l,n,a;m,p,@ transforms as a second-rank tensor. How- 
ever, for a symmetry operation, the transformed matrix 
must be the same as the (unrotated) matrix of force con- 
stants between the pair of ions which are mapped to (I, n) 
and (m,p). Hence, new elements of the matrix of force 
constants can be obtained by application of these point 
symmetries. Translational symmetries can be identified 
and exploited in a similar fashion. 

The matrix of force constants should be symmetric fjlf l 
and Newton's third law must be satisfied: if an ion is dis- 
placed slightly then the restoring force on that ion must 
be equal and opposite to the total force on all of the other 
ions. Hence we must have: 

$l,n,a;m : p,{3 — 0m,p,/3;/,n,a (1) 
4>l,n,a;l,n,P = — J]] 4>l,n,a;m,p,0 ■ (2) 

The force constants are obtained from separate numer- 
ical calculations; hence small violations of these require- 



ments may occur. These two conditions are therefore al- 
ternately imposed on the matrix of force constants until 
further application leaves the matrix unchanged 11 . 

The next stage of the calculation involves the construc- 
tion of dynamical matrices [ p"3[ for various wavevectors 
in the first Brillouin zone. Diagonalization of the dynam- 
ical matrix for a given wavevector gives the spectrum of 
corresponding eigenfrequencies. Strictly, these are only 
exact when the wavelengths are commensurate with the 
dimensions of the supercell pi]] . However, provided that 
the resulting dispersion curves are smooth, it may be as- 
sumed that the interpolation errors are negligible. 

Cochran and Cowley |l4|] have shown that the elements 
of the dynamical matrix for an ionic crystal can be writ- 
ten as the sum of a term that behaves analytically as the 
wavevector tends to zero and a term that is non-analytic 
at the zone center. The latter term vanishes as the 
boundary of the Brillouin zone is approached. This term 
arises because longitudinal optic (LO) phonons cause an 
electric polarization field to be set up within the crystal 
as the oppositely-charged ions are displaced in opposite 
directions. The resulting long-range interactions cannot 
be calculated within the framework we have described so 
far because of the limited size of the simulation supercell. 
At the zone center itself, the LO phonon sets up a uni- 
form electric polarization that is incompatible with the 
periodic boundary conditions on the supercell. 

Cochran and Cowley's expression for the dynamical 
matrix is: 



Aire 1 



Q|k| 2 /MJWp 

3 \ * 

^ ^ ^7-^n,7,a(k) J 6q (k) 



Vy=l 

X ^fc 7 Z p , 7 ,g(k) 
Vy=l J 



(3) 



where k is the wavevector, e is the electronic charge, f2 
is the volume of the unit cell, M n is the mass of ion 
n , Z n ^ aj i3(\i) is the Born effective charge tensor for ion 
n and eo(k) is the electronic (frequency dependent) di- 
electric function. The first term on the right-hand side 
is the component of the dynamical matrix that is ana- 
lytic as k — > 0, while the second term is the non-analytic 
part due to macroscopic polarization effects. We use our 
matrix of force constants evaluated using the Hellmann- 
Feynman theorem in a cubic supercell to evaluate the 
analytic part as: 
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where R m is the position vector of unit cell m. 

Following Parlinski et al |L2| we assume that the sec- 
ond term on the right-hand side of Equation || falls off 
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from its value at the Brillouin zone center with a Gaus- 
sian profile. For wavevectors in the first Brillouin zone 
this term is: 

47re 2 ( ^ , \* 

x ^fc 7 Z p , %/3 (0)^ xe^(^) , (5) 

where k 1 / 2 is the distance from the center to the bound- 
ary of the Brillouin zone along the k x -, k y - and /in- 
directions and po is a parameter determining the rate 
at which the term falls off as the edge of the Brillouin 
zone is approached. Following Parlinski, we set po = 1.2. 

The frequency density-of-states function is evaluated 
using the method of Swift in which the Brillouin 
zone is sampled using Monte Carlo methods. For a sin- 
gle harmonic mode of frequency to, the Hclmholtz free 
energy is given by: 
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where h is the Dirac constant, fc^ is Boltzmann's con- 
stant, T is the temperature and [3 = 1/ksT. Hence, 
by numerically integrating the product of the specific 
density-of-states with the mean free energy of a normal 
mode, the specific lattice thermal free energy can be cal- 
culated for a range of temperatures. The dependence on 
density is found by interpolating between the results at 
different lattice parameters. 



B. Polymorphism 

C. The Gibbs free energy 

We have calculated the Helmholtz free energy f(v,T) 
as a function of temperature T and specific volume v 
(the reciprocal of the density). However, the appropri- 
ate thermodynamic potential for constructing the (p, T)- 
phase diagram and evaluating the polymorphic equation 
of state is the specific Gibbs free energy g(p, T), where p 
is the pressure. The Gibbs free energy function for each 
phase can be evaluated using the Legendre transforma- 
tion: 

g(p,T) = f+pv = f-(^\ V. (7) 



D. The phase diagram 

Under conditions of fixed pressure and temperature, 
the system consists entirely of the available phase with 
the lowest Gibbs free energy. Thus the phase diagram in 
(p, T)-space can be evaluated. 



E. Combining phases 



For each pressure and temperature, we may evaluate 
the polymorphic Gibbs free energy g po i y (p, T) as the low- 
est of the Gibbs free energies for each phase. Given this, 
we may carry out a Legendre transformation to the poly- 
morphic Hclmholtz free energy: 
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Differentiating this, we obtain the pressure as a func- 
tion of specific volume and temperature: the desired 
polymorphic equation of state. 



III. EXTENSION OF THE QUASIHARMONIC 
METHOD TO UNSTABLE PHONONS 



A. Analytic model of soft-mode phonons 

In minerals such as perovskites it is possible to de- 
scribe the transition from a high-temperature phase to a 
low-tempcrature phase of lesser symmetry as the "freez- 
ing in" of a finite amplitude of an unstable phonon of 
the high-symmetry phase, plus a finite strain on the unit 
cell. We consider the application of quasiharmonic ideas 
to these materials. 

The simple harmonic model gives a negative energy 
and divergent free energy arising from the unstable 
modes. In reality, the soft-mode phonon is best described 
by a potential double-well with a local maximum at 
the mean structure, corresponding to the high-symmetry 
phase. 

Let Xi be a coordinate describing the structural feature 
involved in the phase transition at a particular wavevec- 
tor. The corresponding normal mode can be modeled by 
considering the dynamics of the set of {xi} moving in 
fixed local potential double- wells |16|. In the harmonic 
limit, normal modes are uncoupled. However, because 
we are considering finite displacements there will in gen- 
eral be coupling between our double-well oscillators, this 
being most pronounced around the phase transition and 
at high temperatures. Coupling can be approximately 
treated by renormalization |16|. 

Much work has been concentrated on the Landau 
model in which the double-well is a free energy in the 
form of a quartic polynomial V(x) — Ax A — B(T)x 2 , 
where B(T) changes sign with temperature through cou- 
pling to other modes. 

Such a polynomial expansion of the total energy is also 
possible, perhaps incorporating still higher-order terms 
P|. However, analytic terms beyond second order imply 
phonon coupling. This is inconsistent with the harmonic 
approximation used to describe non-soft modes: even at 
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high phonon number the normal modes are assumed to 
be harmonic and therefore independent of each other. 

We propose instead to describe the entire soft-phonon 
branch via a double-well of form: 
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V{x) = -moj 2 x 2 + e(e- x /2a ~ - 1). 



(9) 



where e, too and a are wavevector dependent. Provided 
that e > muj^a 2 , there are minima at: 

x = x± = ±J2a 2 log{e/mtu 2 a 2 ), (10) 
separated by a barrier of height: 

AV = V(0) - V(x±) 

= e - muj 2 a 2 (l + log (e/moj 2 a 2 )) . (11) 

This form of potential has the advantage of being ap- 
proximately quadratic in both the low-energy and high- 
energy limits. Specifically, for the low-energy case at 
x = x±, we have: 



d 2 V 
dx 2 
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which is equivalent to an harmonic oscillator of fre- 
quency: 



j'o = ^2uj 2 \og(e/muj 2 a 2 ). 
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On the other hand, for the high-energy case, the potential 
approximates that of an harmonic oscillator of frequency 

LOq. 

Soft modes do not usually show abnormal dependence 
on temperature — except in the vicinity of the phase tran- 
sition. Therefore we expect our model to be more widely 
applicable than models where soft modes are treated as 
quartic and other modes as harmonic. 

The (imaginary) harmonic frequency uj c about the (un- 
stable) center of the well is given by: 
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where H(p,x) = p 2 /2m + V(x) is the Hamiltonian of 
the isolated mode as a function of x and p, the canoni- 
cal momentum conjugate to x. (3 — l/kgT where ks is 
Boltzmann's constant and T is the temperature. 

For a given energy E, the frequency of our isolated 
mode can be evaluated using the action-angle method. 
The action variable is: 
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pdx. 
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If the mode has energy E > e, so that there is sufficient 
energy to cross the barrier each libration, we find that: 



2fj f x M{E)/a I ■ 

j = — J y 2mE ~ m 2 Lola 2 z 2 - 2mee~ z2 / 2 dz, 
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where xm is the positive solution of V(x) — E. On the 
other hand, if E < e, so that the motion is confined to 
one side of the double- well, we find that: 
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(18) 



where xm is the greater of the two positive solutions and 
x m is the lesser. 

In either case, the corresponding frequency u> may be 
evaluated using Equation n9[ 



B. Isolated double- well oscillators 

We now consider the problem of motion in an isolated 
potential double-well. 

1. Classical solution 

To evaluate the mechanical energy, we assume the 
mode is in thermal contact with a heat bath at the ap- 
propriate temperature. The mean energy is given by: 



dE' 



(19) 



Taken together, the results of Equations |T^, [l?], [T^ and 
|l9l allow us to calculate numerically the frequency of an 
isolated oscillator moving with the mean thermal energy 
as a function of temperature. Example results are shown 
in Figure [l[ Typical "soft-mode" behavior is observed, 
with the frequency dropping to zero in a cusp at the 
transition temperature. This simple approach was used 
in early studies of MgSi0 3 @. 
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2. Quantum solution 

The energy eigenfunctions of a particle moving in a 
symmetric potential must be either symmetric or anti- 
symmetric. Furthermore, by considering building up the 
Gaussian barrier adiabatically, it is clear that the sym- 
metry of the eigenfunctions must be the same as for those 
of the harmonic oscillator. 

The definite symmetry of the wavefunctions leads to a 
"paradox" for wells of finite separation. If we know the 
energy of the oscillator then it is in an energy eigenstate 
and the wavefunction is either symmetric or antisym- 
metric. Hence the probability distribution is symmetric 
about the center of the double- well and we cannot mean- 
ingfully say which side the oscillator is confined in, even 
if its energy is much less than the barrier height. Thus it 
is not conceptually clear that equating the mean thermal 
energy with the barrier height gives the correct transition 
temperature. We discuss this further in Section [II D| . 

The Hamiltonian operator for a particle moving in the 
quadratic potential without the additional Gaussian po- 
tential is: 



2m dx 2 
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The well-known energy eigenvalues and eigenfunctions 
for the time-independent Schrodinger equation H°(f> n = 
E°<j> n are: 



E° n = (n + l/2)ftw 



and 



(21) 



(22) 



for n G Mi, where H n {x) is the nth Hermite polynomial. 

The Hamiltonian operator for a particle moving in the 
double-well is H = H° + Vi where V x = eie'* 2 / 2 * 2 - 1) 
is the extra Gaussian term. Let the eigenfunctions and 
eigenenergies of the full Hamiltonian be rp n and E n . 

The eigenfunctions of the simple harmonic oscillator 
(SHO) are chosen as the basis of wavefunction space. 
This choice makes the computation particularly simple, 
as will be seen below. 

The matrix elements of the Hamiltonian with respect 
to our chosen basis are: 

(&|ff&) = <&l#Vj) + <&l*Wj) 

= (i + l/2)huj Si } j 
e 



2(' 1 +J')/ 2 V^ 



e~ Kz H i (z)H J (z)dz, 
(23) 



where K = h/2muioa 2 + 1. Note that the matrix is real 
and symmetric. 

The Hermite polynomials satisfy Hi{—x) — 
{-lfHi(x), so that (falVifa) = if i + j is odd. If 



i + j is even then the integrand is an even function. 
Hence, in this case: 
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e- Kz 'H i {z)H J (z)dz. 



(24) 



The eigenvalues of the matrix of the Hamiltonian are 
the allowed energy levels. For energies that are large 
compared with the barrier height the particle will spend 
most of its time away from the center of the well; hence 
we expect that the system will behave as a SHO in this 
limit. Indeed, it is clear from Equation |24] that the ele- 
ments of the matrix (4>i\Vi(f>j) fall off rapidly as i and j 
increase. Hence, for large i or j, the eigenvalues of the 
Hamiltonian tend to to those of the SHO. Thus we only 
need to diagonalize the upper left-hand corner (say, the 
(n c + 1) x (n c + 1) submatrix) in order to obtain the first 
n c + 1 energy levels Eq to E Uc . Beyond n c the eigen- 
values may be taken to be those of the SHO. The com- 
parative ease with which the Hamiltonian matrix can be 
diagonalized is one of the advantages of the quadratic- 
plus-gaussian double-well over the quartic double-well 
potential, although there is no analytic form equivalent 
to equation ^ 

Assuming that n c is sufficiently large, the canonical 
partition function can be written as: 



n=0 
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(25) 



The free energy of the double-well oscillator can then 
be evaluated as: 



F 1 = -k B Tlog(Z). 



(26) 



C. Practical implementation in the quasiharmonic 
method 



Having proposed that each soft mode at a given 
wavevector be described by a double-well of the form 
given in Equation ^, we now describe how the param- 
eters e, a and loq can be determined. Note that if the 
{xi} are mass- reduced phonon coefficients then we may, 
without loss of generality, set m = 1. 

Consider the phonon dispersion curve of a crystal 
structure in which imaginary frequencies are present. 
Those branches that remain real throughout the whole 
of the Brillouin zone are treated as harmonic and, for 
each mode at each wavevector, Equation ^ may be used 
to find the corresponding free energy. For those branches 
that are imaginary in some region of the Brillouin zone, 
however, we propose the following treatment: 
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1. At each symmetry point of the Brillouin zone the 
eigenvector corresponding to the relevant mode 
should be evaluated and the displacement pat- 
tern frozen into the crystal. Ab initio techniques 
can then be used to find the corresponding low- 
symmetry structure, if desired. Using three total 
energy calculations with different amplitudes of the 
soft phonon frozen into the structure 1 , we may eval- 
uate the parameters e, a and wo of the double- well. 
Note that it is possible to fit equation || to every 
branch even if the mode is not imaginary since that 
Equation ^ does not necessarily describe a double- 
well. In practice the harmonic approximation is 
used for all-real branches, it is equivalent to setting 
e = 0. 

2. For each branch we use our results for the double- 
well parameters at the symmetry points to con- 
struct interpolating polynomials over the whole of 
the Brillouin zone for the e and a parameters. 

3. For any wavevector in the Brillouin zone, we may 
find the spectrum of corresponding (possibly imag- 
inary) frequencies. Provided we know to which 
branch these modes belong, we have sufficient in- 
formation to determine the parameters of the ap- 
propriate double-well for each mode, e and a are 
found by interpolating to our wavevector and the 
unstable frequency gives oj c (see Equation |l^) , from 
which we may find u> 2 — u 2 + e/ma 2 . 

4. Hence, for any given wavevector, the free energy of 
each mode, whether harmonic or soft, can be eval- 
uated. These free energies can be summed to give 
the free energy contribution from all modes at the 
given wavevector. 

5. The free energy can then be integrated over all 
wavevectors in the Brillouin zone to give the total 
lattice thermal free energy. By using a grid-based 
scheme to integrate over an irreducible wedge of 
the zone and by making use of the continuity of 
the gradient of each branch, it is possible to keep 
track of which branch is which — necessary if the ap- 
propriate values of e and a are to be interpolated in 
the pressence of imaginary branch crossings. The 
problem of interpolation over an irreducible wedge 
of the Brillouin zone has been studied extensively 
in the context of electronic eigenvalues: see, for ex- 
ample, Reference Q. 

The high-symmetry dynamically stabilized phase and 
the low-symmetry "frozen-phonon" phase can now be 



trea ted a s being distinct. Hence the methodology of Sec- 
tion 



II B can be applied to find the phase diagram. 



D. Interpretation of the soft mode transition 

Consider an isolated symmetric double- well oscillator. 
For the probability density to be asymmetric — necessary 
if we are to meaningfully say that the particle is in one 
well or the other — we must have a mixture of symmetric 
and antisymmetric energy eigenstates. Therefore we can- 
not simultaneously know the energy of our particle and 
which well it is in unless we break the symmetry (e.g. 
by allowing the crystal to distort under phonon-strain 
coupling). 

If the oscillating particle's wavefunction is a superpo- 
sition of different energy eigenfunctions then the expan- 
sion coefficients will evolve in time (provided the sys- 
tem remains both undisturbed and unobserved) accord- 
ing to the Schrodinger equation. Hence the quantum 
mechanical expectation value of the particle's position, 
(x) changes in time. 

The time-dependent wavefunction can be written as: 
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1p n (x). 
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Hence the expectation of x can be written as: 
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where uj nrn = \E n - E m \/% and r\ nm = (arg(c m ) - 
arg(c n ))sgn(w„ m ). Note that we use the fact that 
(ipnlx^m) = if tp n and ip m are either both odd or 
both even since x itself is odd. We also use the fact 
that {ip n \xip m ) = (ip m \xip n ). 

We assume that the energy levels are initially pop- 
ulated according to Boltzmann statistics; thus |c„' 2 
z -i e -/3B„ 5 where p = l j kBT and z = J2™ =0 e~ 

the canonical partition function. 
So we have: 
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(x) = T nm cos (a; 



(29) 



n even m odd 



1 In fact, if we have calculated the total energy of the high- 
symmetry phase (corresponding to zero amplitude of the 
phonon) then we only need to carry out a further two frozen 
phonon calculations. 
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where T nm = 2Z- 1 e~ ti{ - E ^ +E ^/ 2 {%l> n \x%l; m ) is the ampli- 
tude of the sinusoidal component in the expansion of (x) 
with frequency u> nm . 

For the harmonic oscillator potential, the frequen- 
cies of the oscillations in (x) are of the form u) nm — 
\E n — E m \/h = \n — m|u>o. Thus the lowest oscillation 
frequency is luq. For the symmetric double- well, how- 
ever, we end up with a set of pairs of energy levels that 
are very close to each other (becoming degenerate in the 
limit that the barrier height goes to infinity) . These give 
rise to oscillation frequencies very much lower than u . 

As the temperature is increased, higher frequency com- 
ponents have larger T coefficients. We suggest that the 
soft mode phase transition be judged to occur when the 
frequency with the highest coefficient exceeds the fre- 
quency of the experimental probe. When this has hap- 
pened, the predominant sinusoidal component of (x) has 
a frequency higher than can be measured by the exper- 
imental probe, and so it appears to the experimenter 
that (x) = 0. Below this temperature, measurements 
of (x) will tend to find it in one well or the other. This 
defin ition is different from the polymorphic one (Section 
[IB) because of the contribution to the free energy from 
the symmetry-breaking distortion of the lattice that in- 
evitably accompanies the transition. In particular, the 
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E = 0.l eV 
e = 0.2 eV 
E = 0.2972 eV 
E = 0.4 eV 



1. The free energy of the soft-mode phase, calculated 
as above, expanded about an unstable frozen-ion 
structure without strain-phonon coupling. 

with 

2. The free energy of the low symmetry phase, cal- 
culated by expanding about the minimum of total 
energy. 

Typically, the former will have higher entropy (sam- 
pling from both wells) while the latter has lower energy. 



E. Absorption of low-frequency photons 

Figure ^| shows the energy difference between neigh- 
boring energy levels (E n — E n -i)/fi plotted against the 
mean of the two energies (E n + E n -i)/2 for the quan- 
tum double-well oscillator. Absorption of photons at fre- 
quency (E n — E n -i)/h is symmetry-allowed. 

For energies in excess of the barrier height the fre- 
quency (E n — E n ^i)/h is virtually identical to the classi- 
cal frequency for energ y (E n + E n -i)/2, obtained using 



the method of Section III B 1 . In the very high-energy 



limit the frequency behaves as that of the quadratic po- 
tential well without the Gaussian barrier. 

For energies less than the barrier height the energy 
levels tend to degenerate pairs of levels. The frequencies 
given by the difference between the energy levels of neigh- 
boring pairs again correspond to the classical frequencies. 
However, the pairs of almost-degenerate eigenstates im- 
ply the existence of very low-frequency absorption peaks. 
(These are the very low frequencies that alternate with 
the classical frequencies below the transition energy in 
Figure ||.) It should be noted that these frequencies 
are not associated with the normal modes of the low- 
symmetry phase and do not, therefore, contribute to the 
quasiharmonic thermal energy. They are a feature of the 



1000 
Temperature / K 



FIG. 1. Classical frequency of a double- well oscillator in 
thermal contact with a heat bath plotted against tempera- 
ture for various barrier heights e. The other double-well pa- 
rameters are: m = 1, ujq = 0.0691 eV 1 ^ 2 A -1 amu -1 ' 2 and 
a = 1.866 amu 1 ' 2 A. The frequency falls to zero in a cusp at 
the phase transition. At e = 0.2972 eV, the double-well pa- 
rameters are appropriate for the double-well describing the 
orthorhombic-tetragonal transition in MgSiOs at zero pres- 
sure [^|. Thus this model predicts a soft mode transition 
temperature of 2609 K, if the coupling of the soft phonon to 
strain is neglected. 

Thus the first-order transition is determined by com- 
paring: 



- Classical frequency 
Difference between neidihrniiiL' encr^\ lc\ Ji \ kled hy Dir; 




Oscillator energy / cV 

FIG. 2. Thick line: Frequency of classical isolated dou- 
ble-well oscillator against energy. The double-well parame- 
ters are as for Figure [j] with e = 0.2972 eV. Fine line: Energy 
differences (Ei — Ei-i)/h against energy (Ei + Ei-i)/2h for 
quantum oscillator. 
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IV. APPLICATION OF QUASIHARMONIC 
METHODS TO PERICLASE 

A. Computational details 

1. The cold-curve 

For each lattice parameter the total energy is evaluated 
using the CASTEP software package ]l9[ , which utilizes 
density-functional theory in the generalized gradient ap- 
proximation (GGA) (20). The ionic cores are accounted 
for using ultrasoft pseudopotentials |2l[| in Kleinman- 
Bylander form [Q. The wavefunctions of the valence 
electrons are expanded in a plane wave basis set up to an 
energy cutoff of 540 eV. 

For the Bl phase, the simulation cell consists of a sin- 
gle cubic unit cell. The Brillouin zone is sampled at 20 
special points generated from an 8 x 8 x 8 mesh using 
the Monkhorst-Pack scheme E3|. For the B2 phase, the 
simulation cell is a single cubic primitive cell. The Bril- 
louin zone is sampled at 35 special points from a 9 x 9 x 9 
mesh. In each case the point symmetries of the crystal 
are enforced p4|| . 

The equilibrium lattice parameter of the Bl phase at 
zero external pressure (which corresponds to the min- 
imum of the cold-curve) calculated using CASTEP is 
a = 4.259 A, which may be compared with an exper- 
imentally determined parameter a — 4.2115(1) A p8[ . 
The difference between the theoretical and experimental 
values is about 1%. 



2. Determination of the matrix of force constants 

The supercells simulated to determine the matrices of 
force constants for the Bl phase consist of 2 x 2 x 2 cubic 
unit cells (64 atoms). Thus the interactions between a 
given ion and its third-closest shell of neighbors are in- 
cluded in our calculations. For the B2 phase, the cells 
used consist of 2 x 2 x 2 cubic primitive unit cells (16 
atoms). In these supercells the crystal symmetry is such 
that only two ionic displacements are required to com- 
plete the entire matrix of force constants. The plane- 
wave cutoff energy is 540 eV and the Brillouin zone is 
sampled at 6 special points from a 4 x 4 x 4 mesh. In each 
simulation the ion displaced from equilibrium is moved 
by 0.4% of the lattice parameter. 

As demonstrated by Parlinski |Q it is possible to cal- 
culate the Born effective charge tensors from first princi- 
ples using simulations of elongated supercells. Note that 
because of the symmetry of the Bl and B2 phases, the 
Born effective charge tensors are isotropic (so Z n ^ a ^ = 
Z n 5 a ,f3)- Furthermore, the sum of the Born effective 
charges over the ions in a unit cell must be zero 
(so Zjvig = —Zo). Hence there is effectively only 
one undetermined parameter in the non-analytic term: 

z Mg (o)/v^(o). 



We choose Z Mg (0) / \/ e (0) = 4.4 to give the LO 
branch in the dispersion curve of Figure |^. Our values 
for the effective charge tensors of the Bl phase are such 
that the calculated LO branch for lattice parameter 4.2 A 
are in reasonable agreement the experimental results of 
Peckham |p5fl . We neglect the variation of the effective 
charge tensors with lattice parameter. 




Waveveclor 



FIG. 3. Dispersion curves for Bl MgO at lattice parame- 
ter 4.2 A. The symmetry points in the dispersion curve are 
(from left to right): V [000], X [001], X [011], Y [000] and 
L [ii §]. Note the LO-TO splitting (At T, LO frequency is 
22.0 THz whereas TO frequency is 12.4 THz) arising from the 
non-analytic term of Equation |E| Also shown are Peckham's 
experimental results for Bl MgO at lattice parameter 4.212 A. 
Note that the reciprocal lattice vectors referred to in the re- 
sults for the Bl phase are those of the cubic unit cell rather 
than the true reciprocal lattice vectors. 



B. Approximations and errors 



1. Errors in ab initio total energy calculations 



Total energy differences between structures calculated 
using density-functional theory in the generalized gradi- 
ent approximation are thought to be reliable to within a 
few percent ]2^| . 

The cutoff energy of the plane- wave basis at 540 eV 
is sufficient for convergence of the total energies of the 
crystals to within 10 /ieV per ion, several orders of mag- 
nitude less than the likely error due to the use of the 
GGA. The Hellmann-Feynman forces are converged to 
within lmeVA -1 , at least two orders of magnitude less 
than the dominant forces arising when an ion is displaced. 
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2. The harmonic approximation 

We investigate the range of validity of the harmonic 
approximation. This is done for the Bl phase with lat- 
tice parameter a = 4.2 A. 

We evaluate the force constant of the restoring force 
on a magnesium ion as it is displaced in the ^-direction. 
The results are shown in Figure [|. It can be seen that the 
force constant starts to increase when the ionic displace- 
ment reaches a max ~ 0.084 A, about 2% of the lattice 
parameter, at which point the potential energy is about 
0.4 eV. Other displacements are similar; hence it is rea- 
sonable to assume that the forces remain linear (and the 
quasiharmonic assumption is valid) for temperatures up 
to several thousand Kelvin. 



Poor convergence of the Hellmann-Feynman forces can 
result in the violation of Newton's third law for the ma- 
trix of force constants. Typically this results in the 
acoustic branches of the dispersion curve failing to pass 
through zero at the center of the Brillouin zone [jn]]. As 
discussed in Section II A 2 , Newton's third law is imposed 
on the matrix of force constants. However, even without 
this, the calculated acoustic branches pass very close to 
zero at the zone center. 

The method by which long-range polarization effects 
are accounted for is also approximate. The effects of this 
are discussed below. 



C. Ab initio phonons 



1. The Bl phase 
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FIG. 4. Graph of restoring force divided by the finite dis- 
placement of a Mg ion in the [100]-direction against that dis- 
placement. Results are appropriate for a 2 x 1 x 1 supercell 
of Bl phase at lattice parameter 4.2 A. The force constant is 
independent of the magnitude of the displacement until the 
displacement is at least 0.084 A. Thereafter the restoring force 
increases faster than linearly with the displacement; the har- 
monic approximation has broken down. 



3. Other sources of error 

Another potential source of error is the limited size of 
the simulation supercells. However, the results shown in 
Table § for the Bl phase make it clear that interactions 
beyond the third-closest shell of neighbors may be safely 
neglected p6| . 

Contributions to the free energy from the thermal ex- 
citation of the valence electrons, from coupled electron- 
phonon excitations and from the equilibrium population 
of defects are thought to be negligible in comparison with 
the frozen- ion and lattice thermal energies ]15[ . 



Inelastic neutron scattering experiments were carried 
out by Peckham |2j| and used to generate dispersion 
curves for the Bl phase p7y . We compare our theoretical 
dispersion curve with these results in Figure p[ (Note 
that our dispersion curve was generated for lattice pa- 
rameter 4.2 A, whereas Peckham's results were obtained 
under ambient conditions where the lattice parameter is 
4.212 A.) Our theoretical results are in reasonable agree- 
ment with experiment. (The lattice parameter of 4.2 A 
corresponds to a pressure of about 7 GPa at zero temper- 
ature.) 

The specific frequency density-of-states function is 
shown in Figure Without the addition of the non- 
analytic term to the dynamical matrix, the longitudi- 
nal optic branch is degenerate with the transverse optic 
branch at the T-point, and this is also shown. Although 
only the LO branch is altered substantially, it can be 
seen that the inclusion of the non-analytic term has a 
significant effect on the density of states. 

We compare sound velocities calculated from our dis- 
persion curves with the experimental results of Reich- 
mann et al ]2§| ] obtained using ultrasonic interferometry. 
Reichmann obtains a P-wave sound speed of 9119 ms -1 
in the [100]-direction whereas our longitudinal-acoustic 
mode has (o!u;LA/dfc[ioo]) k _ = 11367 ms" 1 . 

On the other hand, in the [11 Indirection, the experi- 
mental P-wave velocity is 10125 ms -1 which may be com- 
pared with our theoretical value of 10818 ms^ 1 . 

Thus Reichmann's experimental results show a higher 
degree of anisotropy than do our theoretical results. The 
discrepancy would appear to be (at least partly) caused 
by the imposition of symmetry and Newton's third law 
on the matrix of force constants jpL X[ : if this procedure 
is not carried out then we find that (dwLA/dfc[ioo]) k=0 = 
10770ms- 1 and («iwLA/<%ii]) k=0 = 12651ms" 1 . Al- 
though Reichmann's P-wave velocities are still somewhat 
less than these theoretical velocities, both are now in 
agreement that the velocity in the [100]-direction is less 
than the velocity in the [11 Indirection. 
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With LO phonon correction 

Without LO phonon correction 




f 5 - 




FIG. 5. Specific frequency density-of-states for Bl phase at 
lattice parameter 4.2 A with and without the inclusion of the 
non-analytic term in the dynamical matrix. 



2. The B2 phase 



FIG. 7. Dispersion curves for B2 MgO at a lattice constant 
of 2.7 A. The symmetry points are as for Figure ^j. Note the 
branches of imaginary phonon frequencies, indicating that the 
structure is mechanically unstable at this volume. 



Typical dispersion curves for the B2 phase at lattice 
parameters 2.0 A and 2.7 A are shown in Figures |and|. 
At zero temperature these lattice parameters correspond 
to pressures of 653 GPa and —6 GPa respectively. 

Note the presence of unstable modes at low pressures 
in the dispersion curve of the B2 phase (Figure 0). We 
find that the B2 phase is structurally unstable for pres- 
sures below about 82 GPa. In our calculations for the 
phase boundary in periclase we do not require our novel 
method for dealing with soft modes because the imagi- 
nary frequencies are only found in the B2 phase at pres- 
sures for which the Bl phase is clearly favored. 



D. Ab initio equation of state 

We plot the pressure against specific volume for a range 
of temperatures in Figure ||. This is the desired thermo- 
dynamic equation of state for periclase. Also shown is a 
third-order Birch-Murnaghan equation of state generated 
from the isothermal bulk modulus and its first derivative 
with respect to volume, which were obtained by means 
of ultrasonic sound velocity measurement [ ^9| . 



Temperature OK (Jackson) 




Wavevector 



FIG. 6. Dispersion curves for B2 MgO at a lattice param- 
eter 2.0 A. The symmetry points in the dispersion curve are 
(from left to right): X [§00], V [000], M [§±0], R [§f§], V 
[000] and X [i 00]. 



FIG. 8. Equation of state of periclase: pressure against 
specific volume at various temperatures. The position of the 
phase transition is clearly visible as a kink in the curves. A 
Birch-Murnaghan equation of state generated from Jackson's 
experimental results is also shown. 
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As a test of the validity of our results, the isothermal 
bulk modulus at zero temperature and external pressure 
(where periclase is entirely in the Bl phase) can be com- 
pared with experimental results. We calculate the bulk 
modulus to be — v (dp/dv) T = 155 GPa, whereas an ex- 
perimentally determined value is 160.3 ± 0.3 GPa p9| . 
The theoretical and experimental results differ by about 
3%. 

We may also compare the pressure derivative of the 
bulk modulus at zero pressure and temperature with ex- 
perimental results. The bulk modulus is found to be al- 
most, but not quite, a linear function of pressure. Fitting 
a straight line to the data from —8.3 GPa to 21.12 GPa we 
find that the the gradient is 4.11, in excellent agreement 
with the experimentally determined value of 4.2 ± 0.2 
p9| . It is found that the pressure derivative of the bulk 
modulus decreases slightly as the pressure is increased. 



E. Ab initio phase diagram 

The theoretical phase diagram of solid periclase is 
shown in Figure |^. At pressures and temperatures below 
and to the left of the phase boundary shown in the dia- 
gram, periclase exists in the Bl phase; above and to the 
right of the boundary it exists in the B2 phase. Duffy ]3C| ] 
has shown experimentally that, at room temperature, the 
Bl phase is stable to pressures of at least 227 GPa. This 
lies well within the Bl region of our theoretical phase 
diagram. It can also be noted that the B2 phase is not 
favored at the pressures at which it is structurally unsta- 
ble. 



Also shown in Figure g is the theoretical phase dia- 
gram obtained by Strachan et al []3l| using molecular 
dynamics simulation. Although qualitatively similar, the 
calculated position and orientation of the B1-B2 phase 
boundary is very different. Karki ||] obtained a B1-B2 
transition pressure of 460 GPa at OK, which also differs 
substantially from that of Strachan. Figure [h] shows the 
Gibbs free energy plotted against pressure for the two 
phases at two different temperatures. The difficulty in 
ascertaining the transition pressures at which the curves 
cross is apparent. This consideration will affect all theo- 
retical calculations of the phase diagram of Dericlasc. 



T^7 




600 



FIG. 10. Gibbs free energy plotted against pressure for the 
Bl and B2 phases of periclase at two different temperatures. 
The transition pressure at each temperature corresponds to 
the point where the curves cross. The similarity of the curves 
for each phase means that small errors in the energies lead to 
large uncertainties in the transition pressures. 
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FIG. 9. Theoretical phase diagram of periclase. The Bl 
phase region is below and to the left of the B1-B2 coexis- 
tence line. Also shown are other theoretical results [p|,[3l|, 
including a theoretical Bl-liquid phase boundary. 



We confirm the difficulty in locating the transition 
by attempting to reproduce Karki's zero-temperature re- 
sults in which zero-point lattice vibrational energy is ne- 
glected. We simply calculate the enthalpy against pres- 
sure for the two phases using CASTEP. The results are 
shown in Figure 0. We find a transition pressure of 
GPa, different from that of Karki (460 GPa). The possi- 
bility of a metallization transition lowering the energy of 
the B2 phase was investigated but found not to occur at 
relevant pressures. 



ff 



Enthalpy against pressure for B 1 and B2 MgO 

Curves cross at 663.975 GPii 




FIG. 11. Enthalpy h — e + pv plotted against pressure 
for the Bl and B2 phases of MgO at temperature K (note 
that zero-point energy is not included in the results shown in 
this graph). The curves cross (and so the phase transition is 
predicted to occur) at 664 GPa. These results are obtained 
directly from the results of CASTEP. Note that, again, the 
two curves are extremely close when they cross; thus the po- 
sition of the theoretical transition point is sensitive to the 
simulation parameters. 

We consider the fractional error in unit cell volume for 
the B2 phase to which the difference between the specific 
enthalpies of the Bl and B2 phases at 451 GPa corre- 
sponds. For given unit cell enthalpy for the B2 phase, we 
find the fractional volume difference as: 



A^ 



V 



L B2 



h B i 



IB2 



(30) 



where h B i = -34.4043 eVA" 3 and h B 2 
-34.3978 eVA" 3 are the specific enthalpies of the two 
phases at 451 GPa. Thus a discrepancy of AV/V x 
10~ 4 « 0.02% in the volume of either phase would cor- 
respond to a 200GPa change in the transition pressure. 
This illustrates the sensitivity of the transition point 
to the details of how the total energy calculations are 
carried out. 

We use the generalized gradient approximation and ul- 
trasoft pseudopotentials p0| , pl| , whereas Karki used the 
local density approximation and Q c tuned pseudopoten- 
tials |3^,|33|]. The difference between calculated volume 
using these two methods is typically of the order of 5%; 
hence this is likely to be responsible for the difference 
between our results and Karki's. 

The Clausius-Clapeyron equation for the coexistence 
line between two phases in (p, T)-space is: 



dp 
dT 



As 
Av' 



(31) 



where p(T) is the coexistence line and As and Av are 
the specific entropy and volume differences between the 
two phases across the line. At zero temperature the 
entropy of the two phases should be zero by the third 
law of thermodynamics which is valid within our method 



(though not, e.g. in classical MD Q); therefore, pro- 
vided the phases have different densities, the coexistence 
line should satisfy dp/dT = at T = 0. Figure ||.) does 
not appear to satisfy this requirement. 

The explanation for this apparent discrepancy lies with 
the fact that the densities of the two phases are very sim- 
ilar (and converging) at their predicted zero-temperature 
transition point: the specific volumes of the Bl and B2 
phases are 0.202 A 3 amu" - 1 and 0.198A 3 amu 1 respec- 
tively. 



V. CONCLUSIONS 

We have described an extension to the quasiharmonic 
method that allows the free energy contribution from 
"soft" phonons in dynamically stabilized crystals to be 
evaluated. Our approach is based on a form of potential 
double-well different to that used in previous work on 
soft phonons: a parabola-plus-Gaussian form that has 
the advantage of being harmonic in both the low- and 
high-temperature limits. 

We argue that the first-order nature of the phase tran- 
sitions found using our extended quasiharmonic method 
arises because of the coupling of the relevant phonon to 
strain in the crystal. Without this coupling, the tran- 
sition would be second-order. We have suggested a cri- 
terion for judging when a second-order soft-mode phase 
transition has occurred, taking into account the quantum 
mechanical nature of the problem. 

At energies less than the height of the central barrier 
in our symmetric potential double-well, the allowed en- 
ergy levels consist of near-degenerate pairs. Hence we 
suggest there must exist extremely low-frequency pho- 
ton absorption peaks for soft-mode materials in their 
low-temperature phase, corresponding to photon-induced 
transitions between such pairs. 

We have evaluated the equation of state and the phase 
boundary for the B1-B2 transition in periclase using ab 
initio calculations in the quasiharmonic approximation. 
We predict that this transition will occur, but that it 
is well outside the ranges encountered inside the Earth. 
Locating the B1-B2 phase boundary with precision is dif- 
ficult, however, because the Gibbs free energy curves of 
the two phases are very similar when they cross. 
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Position of Mg atom 


Magnitude of [001]-component 




of force / eVA -1 


(0,0,0) 


0.13424 


(0.5,0.5,0.0625) 


0.03356 


(0,0,0.125) 


0.00981 


(0.5,0.5,0.1875) 


0.00006 



TABLE I. Magnitude of the component of force in the 
[001]-direction when the Mg ion at (0,0,0) is displaced in the 
[001]-direction by 0.05% of the length of a l/\/2 x l/\/2 x 8 
supercell. Coordinates are given as fractions of the supercell 
dimensions. (These results are for Bl MgO at lattice param- 
eter 4.2 A.) 
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